Set Up¶

In [1]:
# Import Libraries

import json                 # Reading and writing JSON
from glob import glob       # File pattern matching
import os                   # Operating system access
import pathlib              # File navigation
import holoviews as hv      # Interactable plot
import panel as pn          # Plot formatting
pn.extension()              # Plot formatting
hv.extension('bokeh')       # Interactable plot

import earthpy              # Work with geospatial data
import xarray as xr         # Multi-dimensional arrays
import hvplot.xarray        # HVPlot and xarray
import rioxarray as rxr     # Raster for xarray

import matplotlib.pyplot as plt    # Plotting
import pandas as pd         # Work with tabular data
import geopandas as gpd     # Wowk with geospatial data
import hvplot.pandas        # Plot geospatial data

Download Data¶

The dataset is saved as 'Gila River Vegetation.'

In [2]:
project = earthpy.Project(
    'Gila River Vegetation', dirname='gric-vegetation-project')
project.get_data()

The 2020 American Indian Tribal Subdivision National (AITSN) borders from the 2020 US census are included in the data downloaded. More information about the shapefiles can be found here.

In [3]:
# Load in the boundary data
aitsn_gdf = gpd.read_file(project.project_dir / 'tl_2020_us_aitsn')

# Check that it worked
print(type(aitsn_gdf))

aitsn_gdf.head()
<class 'geopandas.geodataframe.GeoDataFrame'>
Out[3]:
AIANNHCE TRSUBCE TRSUBNS GEOID NAME NAMELSAD LSAD CLASSFP MTFCC FUNCSTAT ALAND AWATER INTPTLAT INTPTLON geometry
0 2430 653 02419073 2430653 Red Valley Red Valley Chapter T2 D7 G2300 A 922036695 195247 +36.6294607 -109.0550394 POLYGON ((-109.2827 36.64644, -109.28181 36.65...
1 2430 665 02419077 2430665 Rock Point Rock Point Chapter T2 D7 G2300 A 720360268 88806 +36.6598701 -109.6166836 POLYGON ((-109.85922 36.49859, -109.85521 36.5...
2 2430 675 02419081 2430675 Rough Rock Rough Rock Chapter T2 D7 G2300 A 364475668 216144 +36.3976971 -109.7695183 POLYGON ((-109.93053 36.40672, -109.92923 36.4...
3 2430 325 02418975 2430325 Indian Wells Indian Wells Chapter T2 D7 G2300 A 717835323 133795 +35.3248534 -110.0855000 POLYGON ((-110.24222 35.36327, -110.24215 35.3...
4 2430 355 02418983 2430355 Kayenta Kayenta Chapter T2 D7 G2300 A 1419241065 1982848 +36.6884391 -110.3045616 POLYGON ((-110.56817 36.73489, -110.56603 36.7...

Through some trial and error, I found latitude and longitude boundaries to focus the plot on the Gila River Indian Community. The AIANNHCE value 1310 can then be used to create a new GDF containing just the GRIC Boundaries.

In [4]:
# Generate plot of surrounding area with Tribal locations highlighted
aitsn_gdf.hvplot(
    geo=True, tiles='EsriImagery',
    frame_width=500,
    legend=False, fill_color=None, edge_color='white',
    hover_cols='all',
    xlim=(-112.4, -111.4),
    ylim=(32.8, 33.6)
)
WARNING:param.main: edge_color option not found for polygons plot with bokeh; similar options include: ['muted_color', 'bgcolor', 'line_color']
Out[4]:
In [5]:
# Select and merge the subdivisions we want (1310)
gric_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='1310'].dissolve()
# Plot the results with web tile images
gric_gdf.hvplot()
Out[5]:
In [6]:
# Create the boundary GDF
boundary_gdf = aitsn_gdf.loc[aitsn_gdf.AIANNHCE=='1310'].dissolve()

# Plot the results with web tile images
boundary_gdf.hvplot(
    geo=True, tiles='EsriImagery',
    fill_color='green', line_color='black',
    fill_alpha=.3,  # must be b/w 0-1
    title='Gila River Indian Community',
    frame_width=500)
Out[6]:

Set up NDVI data¶

Create and combine a list of the raster files with the unique date for each raster now included as an extra dimension to the (x,y) spatial data.

In [7]:
# Get a sorted list of NDVI tif file paths
ndvi_paths = sorted(list(project.project_dir.rglob('*NDVI*.tif')))

# Display the first and last three files paths to check the pattern
ndvi_paths[:3], ndvi_paths[-3:]
Out[7]:
([PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001145000000_aid0001.tif'),
  PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001161000000_aid0001.tif'),
  PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001177000000_aid0001.tif')],
 [PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022209000000_aid0001.tif'),
  PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022225000000_aid0001.tif'),
  PosixPath('/workspaces/data/gric-vegetation-project/gila-river-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022241000000_aid0001.tif')])
In [8]:
# Loop through each NDVI image to add date dimension

doy_start = -25
doy_end = -19

ndvi_das = []

for ndvi_path in ndvi_paths:
    # Get date from file name
    doy = ndvi_path.name[doy_start:doy_end]
    date = pd.to_datetime(doy, format='%Y%j')

    # Open dataset
    da = rxr.open_rasterio(ndvi_path, mask_and_scale=True).squeeze()

    # Add date dimension and clean up metadata
    da = da.assign_coords({'date': date})
    da = da.expand_dims({'date': 1})
    da.name = 'NDVI'

    # Prepare for concatenation
    ndvi_das.append(da)

len(ndvi_das)
Out[8]:
154

Combine the list of rasters into a new data xarray "sliced" by the dates.

In [9]:
# Combine NDVI images from all dates
ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])

ndvi_da
/tmp/ipykernel_1220/646390900.py:2: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
/tmp/ipykernel_1220/646390900.py:2: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
Out[9]:
<xarray.Dataset> Size: 48MB
Dimensions:      (date: 154, y: 203, x: 382)
Coordinates:
    band         int64 8B 1
  * x            (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
  * y            (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
    spatial_ref  int64 8B 0
  * date         (date) datetime64[ns] 1kB 2001-01-14 2001-01-16 ... 2022-01-24
Data variables:
    NDVI         (date, y, x) float32 48MB 0.8282 0.6146 ... 0.2146 0.2085
xarray.Dataset
    • date: 154
    • y: 203
    • x: 382
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      -112.3 -112.3 ... -111.5 -111.5
      array([-112.309375, -112.307292, -112.305208, ..., -111.519792, -111.517708,
             -111.515625], shape=(382,))
    • y
      (y)
      float64
      33.39 33.39 33.38 ... 32.97 32.97
      array([33.388542, 33.386458, 33.384375, ..., 32.971875, 32.969792, 32.967708],
            shape=(203,))
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      GeoTransform :
      -112.31041665660531 0.0020833333331466974 0.0 33.38958333034212 0.0 -0.0020833333331466974
      array(0)
    • date
      (date)
      datetime64[ns]
      2001-01-14 ... 2022-01-24
      array(['2001-01-14T00:00:00.000000000', '2001-01-16T00:00:00.000000000',
             '2001-01-17T00:00:00.000000000', '2001-01-19T00:00:00.000000000',
             '2001-01-20T00:00:00.000000000', '2001-01-22T00:00:00.000000000',
             '2001-01-24T00:00:00.000000000', '2002-01-14T00:00:00.000000000',
             '2002-01-16T00:00:00.000000000', '2002-01-17T00:00:00.000000000',
             '2002-01-19T00:00:00.000000000', '2002-01-20T00:00:00.000000000',
             '2002-01-22T00:00:00.000000000', '2002-01-24T00:00:00.000000000',
             '2003-01-14T00:00:00.000000000', '2003-01-16T00:00:00.000000000',
             '2003-01-17T00:00:00.000000000', '2003-01-19T00:00:00.000000000',
             '2003-01-20T00:00:00.000000000', '2003-01-22T00:00:00.000000000',
             '2003-01-24T00:00:00.000000000', '2004-01-14T00:00:00.000000000',
             '2004-01-16T00:00:00.000000000', '2004-01-17T00:00:00.000000000',
             '2004-01-19T00:00:00.000000000', '2004-01-20T00:00:00.000000000',
             '2004-01-22T00:00:00.000000000', '2004-01-24T00:00:00.000000000',
             '2005-01-14T00:00:00.000000000', '2005-01-16T00:00:00.000000000',
             '2005-01-17T00:00:00.000000000', '2005-01-19T00:00:00.000000000',
             '2005-01-20T00:00:00.000000000', '2005-01-22T00:00:00.000000000',
             '2005-01-24T00:00:00.000000000', '2006-01-14T00:00:00.000000000',
             '2006-01-16T00:00:00.000000000', '2006-01-17T00:00:00.000000000',
             '2006-01-19T00:00:00.000000000', '2006-01-20T00:00:00.000000000',
             '2006-01-22T00:00:00.000000000', '2006-01-24T00:00:00.000000000',
             '2007-01-14T00:00:00.000000000', '2007-01-16T00:00:00.000000000',
             '2007-01-17T00:00:00.000000000', '2007-01-19T00:00:00.000000000',
             '2007-01-20T00:00:00.000000000', '2007-01-22T00:00:00.000000000',
             '2007-01-24T00:00:00.000000000', '2008-01-14T00:00:00.000000000',
             '2008-01-16T00:00:00.000000000', '2008-01-17T00:00:00.000000000',
             '2008-01-19T00:00:00.000000000', '2008-01-20T00:00:00.000000000',
             '2008-01-22T00:00:00.000000000', '2008-01-24T00:00:00.000000000',
             '2009-01-14T00:00:00.000000000', '2009-01-16T00:00:00.000000000',
             '2009-01-17T00:00:00.000000000', '2009-01-19T00:00:00.000000000',
             '2009-01-20T00:00:00.000000000', '2009-01-22T00:00:00.000000000',
             '2009-01-24T00:00:00.000000000', '2010-01-14T00:00:00.000000000',
             '2010-01-16T00:00:00.000000000', '2010-01-17T00:00:00.000000000',
             '2010-01-19T00:00:00.000000000', '2010-01-20T00:00:00.000000000',
             '2010-01-22T00:00:00.000000000', '2010-01-24T00:00:00.000000000',
             '2011-01-14T00:00:00.000000000', '2011-01-16T00:00:00.000000000',
             '2011-01-17T00:00:00.000000000', '2011-01-19T00:00:00.000000000',
             '2011-01-20T00:00:00.000000000', '2011-01-22T00:00:00.000000000',
             '2011-01-24T00:00:00.000000000', '2012-01-14T00:00:00.000000000',
             '2012-01-16T00:00:00.000000000', '2012-01-17T00:00:00.000000000',
             '2012-01-19T00:00:00.000000000', '2012-01-20T00:00:00.000000000',
             '2012-01-22T00:00:00.000000000', '2012-01-24T00:00:00.000000000',
             '2013-01-14T00:00:00.000000000', '2013-01-16T00:00:00.000000000',
             '2013-01-17T00:00:00.000000000', '2013-01-19T00:00:00.000000000',
             '2013-01-20T00:00:00.000000000', '2013-01-22T00:00:00.000000000',
             '2013-01-24T00:00:00.000000000', '2014-01-14T00:00:00.000000000',
             '2014-01-16T00:00:00.000000000', '2014-01-17T00:00:00.000000000',
             '2014-01-19T00:00:00.000000000', '2014-01-20T00:00:00.000000000',
             '2014-01-22T00:00:00.000000000', '2014-01-24T00:00:00.000000000',
             '2015-01-14T00:00:00.000000000', '2015-01-16T00:00:00.000000000',
             '2015-01-17T00:00:00.000000000', '2015-01-19T00:00:00.000000000',
             '2015-01-20T00:00:00.000000000', '2015-01-22T00:00:00.000000000',
             '2015-01-24T00:00:00.000000000', '2016-01-14T00:00:00.000000000',
             '2016-01-16T00:00:00.000000000', '2016-01-17T00:00:00.000000000',
             '2016-01-19T00:00:00.000000000', '2016-01-20T00:00:00.000000000',
             '2016-01-22T00:00:00.000000000', '2016-01-24T00:00:00.000000000',
             '2017-01-14T00:00:00.000000000', '2017-01-16T00:00:00.000000000',
             '2017-01-17T00:00:00.000000000', '2017-01-19T00:00:00.000000000',
             '2017-01-20T00:00:00.000000000', '2017-01-22T00:00:00.000000000',
             '2017-01-24T00:00:00.000000000', '2018-01-14T00:00:00.000000000',
             '2018-01-16T00:00:00.000000000', '2018-01-17T00:00:00.000000000',
             '2018-01-19T00:00:00.000000000', '2018-01-20T00:00:00.000000000',
             '2018-01-22T00:00:00.000000000', '2018-01-24T00:00:00.000000000',
             '2019-01-14T00:00:00.000000000', '2019-01-16T00:00:00.000000000',
             '2019-01-17T00:00:00.000000000', '2019-01-19T00:00:00.000000000',
             '2019-01-20T00:00:00.000000000', '2019-01-22T00:00:00.000000000',
             '2019-01-24T00:00:00.000000000', '2020-01-14T00:00:00.000000000',
             '2020-01-16T00:00:00.000000000', '2020-01-17T00:00:00.000000000',
             '2020-01-19T00:00:00.000000000', '2020-01-20T00:00:00.000000000',
             '2020-01-22T00:00:00.000000000', '2020-01-24T00:00:00.000000000',
             '2021-01-14T00:00:00.000000000', '2021-01-16T00:00:00.000000000',
             '2021-01-17T00:00:00.000000000', '2021-01-19T00:00:00.000000000',
             '2021-01-20T00:00:00.000000000', '2021-01-22T00:00:00.000000000',
             '2021-01-24T00:00:00.000000000', '2022-01-14T00:00:00.000000000',
             '2022-01-16T00:00:00.000000000', '2022-01-17T00:00:00.000000000',
             '2022-01-19T00:00:00.000000000', '2022-01-20T00:00:00.000000000',
             '2022-01-22T00:00:00.000000000', '2022-01-24T00:00:00.000000000'],
            dtype='datetime64[ns]')
    • NDVI
      (date, y, x)
      float32
      0.8282 0.6146 ... 0.2146 0.2085
      units :
      NDVI
      AREA_OR_POINT :
      Area
      array([[[0.8282    , 0.6146    , 0.3796    , ..., 0.1542    ,
               0.1542    , 0.1774    ],
              [0.52349997, 0.32029998, 0.46359998, ..., 0.145     ,
               0.145     , 0.145     ],
              [0.52349997, 0.38509998, 0.4567    , ..., 0.0962    ,
               0.058     , 0.0641    ],
              ...,
              [0.13759999, 0.13759999, 0.1283    , ..., 0.1489    ,
               0.1489    , 0.1568    ],
              [0.1299    , 0.1299    , 0.12379999, ..., 0.1831    ,
               0.1831    , 0.23169999],
              [0.1347    , 0.1309    , 0.1329    , ..., 0.15709999,
               0.2384    , 0.2452    ]],
      
             [[0.5146    , 0.42249998, 0.3706    , ..., 0.1656    ,
               0.1656    , 0.1656    ],
              [0.3543    , 0.4138    , 0.3706    , ..., 0.1275    ,
               0.1275    , 0.0968    ],
              [0.48029998, 0.3543    , 0.3687    , ..., 0.075     ,
               0.075     , 0.0968    ],
      ...
              [0.1303    , 0.1303    , 0.1303    , ..., 0.1739    ,
               0.1739    , 0.1859    ],
              [0.13319999, 0.13319999, 0.1303    , ..., 0.19399999,
               0.19399999, 0.2383    ],
              [0.1662    , 0.13319999, 0.1348    , ..., 0.17289999,
               0.2034    , 0.23249999]],
      
             [[0.5902    , 0.3793    , 0.3145    , ..., 0.15359999,
               0.15359999, 0.1564    ],
              [0.6694    , 0.4878    , 0.32909998, ..., 0.159     ,
               0.159     , 0.1331    ],
              [0.3103    , 0.39409998, 0.4527    , ..., 0.154     ,
               0.1026    , 0.0448    ],
              ...,
              [0.1338    , 0.1338    , 0.1192    , ..., 0.174     ,
               0.174     , 0.163     ],
              [0.12889999, 0.12889999, 0.1205    , ..., 0.2146    ,
               0.2146    , 0.2025    ],
              [0.13069999, 0.1157    , 0.12629999, ..., 0.17359999,
               0.2146    , 0.2085    ]]], shape=(154, 203, 382), dtype=float32)
    • x
      PandasIndex
      PandasIndex(Index([-112.30937498993873, -112.30729165660559, -112.30520832327244,
             -112.30312498993929, -112.30104165660615,   -112.298958323273,
             -112.29687498993985,  -112.2947916566067, -112.29270832327356,
             -112.29062498994041,
             ...
             -111.53437499000816, -111.53229165667501, -111.53020832334187,
             -111.52812499000872, -111.52604165667557, -111.52395832334243,
             -111.52187499000928, -111.51979165667613, -111.51770832334299,
             -111.51562499000984],
            dtype='float64', name='x', length=382))
    • y
      PandasIndex
      PandasIndex(Index([ 33.38854166367555,   33.3864583303424,  33.38437499700925,
              33.38229166367611,  33.38020833034296,  33.37812499700981,
             33.376041663676666,  33.37395833034352,  33.37187499701037,
             33.369791663677226,
             ...
             32.986458330378234,  32.98437499704509,  32.98229166371194,
             32.980208330378794,  32.97812499704565,   32.9760416637125,
             32.973958330379354,  32.97187499704621,  32.96979166371306,
             32.967708330379914],
            dtype='float64', name='y', length=203))
    • date
      PandasIndex
      PandasIndex(DatetimeIndex(['2001-01-14', '2001-01-16', '2001-01-17', '2001-01-19',
                     '2001-01-20', '2001-01-22', '2001-01-24', '2002-01-14',
                     '2002-01-16', '2002-01-17',
                     ...
                     '2021-01-20', '2021-01-22', '2021-01-24', '2022-01-14',
                     '2022-01-16', '2022-01-17', '2022-01-19', '2022-01-20',
                     '2022-01-22', '2022-01-24'],
                    dtype='datetime64[ns]', name='date', length=154, freq=None))
In [10]:
# Plot first and last files side by side to visualize progress
first_ndvi = rxr.open_rasterio(ndvi_paths[0], mask_and_scale=True).squeeze()
last_ndvi = rxr.open_rasterio(ndvi_paths[-1], mask_and_scale=True).squeeze()

fig, axes = plt.subplots(1, 2, figsize = (12, 4))

first_ndvi.plot(ax=axes[0], cmap=plt.cm.PiYG, vmin=-1, vmax=1)
axes[0].set_title("Gila River NDVI - 2001")
boundary_gdf.plot(ax=axes[0], edgecolor='black', facecolor='none',
                  linewidth=1)


last_ndvi.plot(ax=axes[1], cmap=plt.cm.PiYG, vmin=-1, vmax=1)
axes[1].set_title("Gila River NDVI - 2022")
boundary_gdf.plot(ax=axes[1], edgecolor='black', facecolor='none',
                  linewidth=1)
Out[10]:
<Axes: title={'center': 'Gila River NDVI - 2022'}, xlabel='x', ylabel='y'>
No description has been provided for this image

Calculate mean values for NDVI from the beginning of the dataset in 2001 to the water rights settlement in 2004 and an xarray of rasters for each year post 2004.

In [11]:
# Compute the mean NDVI values before water rights were restored
ndvi_pre = (ndvi_da
             .sel(date=slice('2001', '2004'))
             .mean('date')
             .NDVI
)

ndvi_pre.plot()
Out[11]:
<matplotlib.collections.QuadMesh at 0x7c4d6b87b850>
No description has been provided for this image
In [12]:
ndvi_post = ndvi_da.sel(date=slice('2005', '2022')).NDVI

ndvi_post_yearly = ndvi_post.groupby('date.year').mean('date')

print (ndvi_post_yearly)

len(ndvi_post_yearly)
<xarray.DataArray 'NDVI' (year: 18, y: 203, x: 382)> Size: 6MB
array([[[0.5380143 , 0.47134283, 0.4491    , ..., 0.17115715,
         0.17115715, 0.175     ],
        [0.53234285, 0.43379998, 0.44574285, ..., 0.14785714,
         0.14785714, 0.16301428],
        [0.50931424, 0.4504714 , 0.40725714, ..., 0.1275    ,
         0.1056857 , 0.07895714],
        ...,
        [0.18635714, 0.18635714, 0.17112856, ..., 0.17614286,
         0.17614286, 0.17889999],
        [0.17275715, 0.17275715, 0.16979997, ..., 0.21141426,
         0.21141426, 0.2217714 ],
        [0.18247142, 0.17312858, 0.18251428, ..., 0.20607142,
         0.26507142, 0.2588    ]],

       [[0.6633143 , 0.5207143 , 0.4350571 , ..., 0.14278571,
         0.14278571, 0.14188571],
        [0.6846143 , 0.5030286 , 0.35744286, ..., 0.11975714,
         0.11975714, 0.12397142],
        [0.5164143 , 0.42664286, 0.3745286 , ..., 0.13028571,
         0.09927142, 0.07134286],
...
        [0.14237143, 0.14237143, 0.13668571, ..., 0.16885713,
         0.16885713, 0.17524286],
        [0.1333857 , 0.1333857 , 0.13745713, ..., 0.21079998,
         0.21079998, 0.21817143],
        [0.15061429, 0.13785715, 0.13837144, ..., 0.16395713,
         0.21935715, 0.2031714 ]],

       [[0.6115714 , 0.45694283, 0.40044284, ..., 0.16761427,
         0.16761427, 0.16799998],
        [0.5378714 , 0.44529995, 0.35509998, ..., 0.15104285,
         0.15104285, 0.15302856],
        [0.33719996, 0.37427142, 0.4122714 , ..., 0.1349    ,
         0.11072857, 0.09608571],
        ...,
        [0.12699999, 0.12699999, 0.12062857, ..., 0.15535715,
         0.15535715, 0.16347143],
        [0.1187857 , 0.1187857 , 0.11717142, ..., 0.19084285,
         0.19084285, 0.22217143],
        [0.13205715, 0.11748571, 0.12167142, ..., 0.17574285,
         0.19661427, 0.21658571]]], shape=(18, 203, 382), dtype=float32)
Coordinates:
    band         int64 8B 1
  * x            (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
  * y            (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
    spatial_ref  int64 8B 0
  * year         (year) int64 144B 2005 2006 2007 2008 ... 2019 2020 2021 2022
Attributes:
    units:          NDVI
    AREA_OR_POINT:  Area
Out[12]:
18
In [13]:
# Calculate difference between mean ndvi before water rights and each year following the settlement

ndvi_diff = ndvi_post_yearly - ndvi_pre
In [14]:
# Check dimensions and coordinates for calculated difference

print(ndvi_diff.dims)
print(ndvi_diff.coords)
('year', 'y', 'x')
Coordinates:
    band         int64 8B 1
  * x            (x) float64 3kB -112.3 -112.3 -112.3 ... -111.5 -111.5 -111.5
  * y            (y) float64 2kB 33.39 33.39 33.38 33.38 ... 32.97 32.97 32.97
    spatial_ref  int64 8B 0
  * year         (year) int64 144B 2005 2006 2007 2008 ... 2019 2020 2021 2022
In [18]:
# Plot the difference

diff_plot = (
    ndvi_diff.hvplot(x='x', y='y', cmap='PiYG', geo=True, groupby='year',
                     clim=(-0.5, 0.5),   
                     title = 'Gila River Indian Community NDVI\n'
                     'Comparing 2001-2004 to Following Years')
    *
    boundary_gdf.hvplot(geo=True, fill_color=None, line_color='black')
)

diff_plot
Out[18]:
BokehModel(combine_events=True, render_bundle={'docs_json': {'9151993b-ad93-44e6-af5e-51ce56d38647': {'version…

Move the slider to the bottom of the plot and set start to 2005.

In [16]:
years = list(diff_plot.kdims[0].values)
years.sort()     # just to be safe

diff_plot_bottom = pn.panel(
    diff_plot,
    widget_location='bottom',
    widgets={
        'year': pn.widgets.DiscreteSlider(
            name='year',
            options=years,
            value=years[0]
        )
    }
)

diff_plot_bottom
Out[16]:
BokehModel(combine_events=True, render_bundle={'docs_json': {'2b12ca0f-b559-4be4-afb7-ff4da0790ccc': {'version…
In [19]:
# Save the new plot

diff_plot_bottom.save('ndvi_diff_plot_new.html', embed=True)
                                               
WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='e0e0ac3a-719d-4a3e-aed1-9ed229c8e554', ...)